#include "nls_2d.h"


void input(struct param* in){
  char tmp[80],line[80];
  FILE *inp, *log;
  int i = 0;
  
  strcpy (tmp, in->runname);
  strcat (tmp, ".inp");
  inp = fopen(tmp,"r");
  if (inp == NULL ) {
    printf("Input file missing\n");
    exit(1);
  }
  strcpy (tmp, in->runname);
  strcat (tmp, ".log");
  log = fopen(tmp,"w");
  printf ("Running ... %s\n", in->runname);
  while (fgets(line, 80, inp )!= NULL ){
    i++;
    sscanf(line, "%s", tmp);
    if (i==2) {
      (in->height) = atof(tmp);
      fprintf(log, "Height = %s\n", tmp);
    }
    if (i==3) {
      (in->width) = atof(tmp);      
      fprintf(log, "Width = %s\n", tmp);
    }
    if (i==4) {    
      (in->Nx) = atoi(tmp);
      fprintf(log, "Number of points on x axis = %s\n", tmp);
    } 
    if (i==5) {    
      (in->Ny) = atoi(tmp);
      fprintf(log, "Number of points on y axis = %s\n", tmp);
    }
    if (i==6) {
      (in->Lx) = atof(tmp);
    }
    if (i==7) {
      (in->Ly) = atof(tmp);
    }
    if (i==8) {    
      (in->cfl) = atof(tmp);
      fprintf(log, "Computational domain is %fx%f. \nCFL = %s\n",2*(in->Lx), 2*(in->Ly), tmp);
    }
  }
  fclose(inp);
  fclose(log);
};
struct cell* setup(struct param* in){
  int Nx = in->Nx, Ny = in->Ny;
  double Lx = in->Lx, Ly = in->Ly;
  int i, j;
  double x, y;
  double dx = 2*Lx/Nx, dy = 2*Ly/Ny;
  double sigma = in->width;
  double a = in->height;
  struct cell* cll;
  struct mesh* msh;

  cll = (struct cell*)malloc(Nx*Ny*sizeof(struct cell));
  msh = (struct mesh*)malloc(Nx*Ny*sizeof(struct mesh));

  // fill data to points in box
  for (j = 0; j < Ny; j++){
    for (i = 0; i < Nx; i++){
      x = -Lx + i*dx;
      y = -Ly + j*dy;
      cll[i+j*Nx].x = x;
      cll[i+j*Nx].y = y;
      cll[i+j*Nx].u = a*exp(-(x*x+y*y)/(sigma*sigma));
      cll[i+j*Nx].v = 0;  
      x = x + 0.5*dx;
      y = y + 0.5*dy;
      msh[i+j*Nx].x = x;
      msh[i+j*Nx].y = y;
      msh[i+j*Nx].u = a*exp(-(x*x+y*y)/(sigma*sigma));
      msh[i+j*Nx].v = 0;   
    }
  }
  // fill pointers to neighbors in box interior
  for (j = 1; j < Ny-1; j++){
    for (i = 1; i < Nx-1; i++ ){ 
      cll[i+j*Nx].mm = &msh[i-1+(j-1)*Nx];
      cll[i+j*Nx].mp = &msh[i-1+j*Nx];
      cll[i+j*Nx].pm = &msh[i+(j-1)*Nx];
      cll[i+j*Nx].pp = &msh[i+j*Nx];

      msh[i+j*Nx].mm = &cll[i+j*Nx];
      msh[i+j*Nx].mp = &cll[i+(j+1)*Nx];
      msh[i+j*Nx].pm = &cll[i+1+j*Nx];
      msh[i+j*Nx].pp = &cll[i+1+(j+1)*Nx];
    }
  }
  // fill pointers to neighbors along y = -Ly and y = Ly
  for (i = 1; i < Nx-1; i++ ){ 
    cll[i].mm = &msh[i-1+(Ny-1)*Nx];
    cll[i].mp = &msh[i-1];
    cll[i].pm = &msh[i+(Ny-1)*Nx];
    cll[i].pp = &msh[i];

    cll[i+(Ny-1)*Nx].mm = &msh[i-1+(Ny-2)*Nx];
    cll[i+(Ny-1)*Nx].mp = &msh[i-1+(Ny-1)*Nx];
    cll[i+(Ny-1)*Nx].pm = &msh[i+(Ny-2)*Nx];
    cll[i+(Ny-1)*Nx].pp = &msh[i+(Ny-1)*Nx];

    msh[i].mm = &cll[i];
    msh[i].mp = &cll[i+Nx];
    msh[i].pm = &cll[i+1];
    msh[i].pp = &cll[i+1+Nx]; 

    msh[i+(Ny-1)*Nx].mm = &cll[i+(Ny-1)*Nx];
    msh[i+(Ny-1)*Nx].mp = &cll[i];
    msh[i+(Ny-1)*Nx].pm = &cll[i+1+(Ny-1)*Nx];
    msh[i+(Ny-1)*Nx].pp = &cll[i+1];
  }
  // fill pointers to neighbors along x = -Lx and x = Lx
  for (j = 1; j < Ny-1; j++ ){ 
    cll[j*Nx].mm = &msh[Nx-1+(j-1)*Nx];
    cll[j*Nx].mp = &msh[Nx-1+j*Nx];
    cll[j*Nx].pm = &msh[(j-1)*Nx];
    cll[j*Nx].pp = &msh[j*Nx];

    cll[Nx-1+j*Nx].mm = &msh[Nx-2+(j-1)*Nx];
    cll[Nx-1+j*Nx].mp = &msh[Nx-2+j*Nx];
    cll[Nx-1+j*Nx].pm = &msh[Nx-1+(j-1)*Nx];
    cll[Nx-1+j*Nx].pp = &msh[Nx-1+j*Nx];

    msh[j*Nx].mm = &cll[j*Nx];
    msh[j*Nx].mp = &cll[(j+1)*Nx];
    msh[j*Nx].pm = &cll[1+j*Nx];
    msh[j*Nx].pp = &cll[1+(j+1)*Nx]; 

    msh[Nx-1+j*Nx].mm = &cll[Nx-1+j*Nx];
    msh[Nx-1+j*Nx].mp = &cll[Nx-1+(j+1)*Nx];
    msh[Nx-1+j*Nx].pm = &cll[j*Nx];
    msh[Nx-1+j*Nx].pp = &cll[(j+1)*Nx];
  }
  // fill pointers at [-Lx,-Ly]
  cll[0].mm = &msh[Nx-1+(Ny-1)*Nx];
  cll[0].mp = &msh[Nx-1];
  cll[0].pm = &msh[(Ny-1)*Nx];
  cll[0].pp = &msh[0];

  msh[0].mm = &cll[0];
  msh[0].mp = &cll[Nx];
  msh[0].pm = &cll[1];
  msh[0].pp = &cll[1+Nx];

  // fill pointers at [Lx,-Ly]
  cll[Nx-1].mm = &msh[Nx-2+(Ny-1)*Nx];
  cll[Nx-1].mp = &msh[Nx-2];
  cll[Nx-1].pm = &msh[Nx-1+(Ny-1)*Nx];
  cll[Nx-1].pp = &msh[Nx-1];

  msh[Nx-1].mm = &cll[Nx-1];
  msh[Nx-1].mp = &cll[Nx-1+Nx];
  msh[Nx-1].pm = &cll[0];
  msh[Nx-1].pp = &cll[Nx];

  //fill pointers at [-Lx,Ly]
  cll[(Ny-1)*Nx].mm = &msh[Nx-1+(Ny-2)*Nx];
  cll[(Ny-1)*Nx].mp = &msh[Nx-1+(Ny-1)*Nx];
  cll[(Ny-1)*Nx].pm = &msh[(Ny-2)*Nx];
  cll[(Ny-1)*Nx].pp = &msh[(Ny-1)*Nx];

  msh[(Ny-1)*Nx].mm = &cll[(Ny-1)*Nx];
  msh[(Ny-1)*Nx].mp = &cll[0];
  msh[(Ny-1)*Nx].pm = &cll[1+(Ny-1)*Nx];
  msh[(Ny-1)*Nx].pp = &cll[1];

  //fill pointers at [Lx,Ly]
  cll[Nx-1+(Ny-1)*Nx].mm = &msh[Nx-2+(Ny-2)*Nx];
  cll[Nx-1+(Ny-1)*Nx].mp = &msh[Nx-2+(Ny-1)*Nx];
  cll[Nx-1+(Ny-1)*Nx].pm = &msh[Nx-1+(Ny-2)*Nx];
  cll[Nx-1+(Ny-1)*Nx].pp = &msh[Nx-1+(Ny-1)*Nx];

  msh[Nx-1+(Ny-1)*Nx].mm = &cll[Nx-1+(Ny-1)*Nx];
  msh[Nx-1+(Ny-1)*Nx].mp = &cll[Nx-1];
  msh[Nx-1+(Ny-1)*Nx].pm = &cll[(Ny-1)*Nx];
  msh[Nx-1+(Ny-1)*Nx].pp = &cll[0];

  return cll;
};

void output(struct cell* cll, struct param* in){
  int i;
  int Nx = in->Nx, Ny = in->Ny;
  struct mesh* msh = cll[0].pp;
  FILE* fh;
  char tmp[80];

  strcpy(tmp, in->runname);
  strcat(tmp, ".dat");
  fh = fopen(tmp, "w");
  for (i = 0; i < Nx*Ny; i++ ){
    fprintf(fh,"%e\t%e\t%e\n",cll[i].x ,cll[i].y, 
	    sqrt(pow(cll[i].u,2)+pow(cll[i].v,2)));
    fprintf(fh,"%e\t%e\t%e\n",msh[i].x ,msh[i].y, 
	    sqrt(pow(msh[i].u,2)+pow(msh[i].v,2)));
  }
};

void fir_half(struct cell* cll){

}
